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QUALITY CONTROL OF METEOROLOGICAL OBSERVATIONS 
Toke Jayachandran and Richard Franke 



Int roduct ion : A major input component to weather prediction 
systems is the observational data, typically geopotential heights 
and winds from weather observing stations around the world. The 
Navy Operational Global Atmospheric Prediction System (NOGAPS) 
produces weather forecasts using a spectral forecast model. 
Preceding the forecast is a multivariate objective analysis, 
followed by a nonlinear normal mode initialization scheme. The 
multivariate objective analysis scheme uses the weather 
observations to produce analyzed fields with minimum statistical 
error. As is to be expected, some of the observations will turn 
out to be erroneous for various reasons such as the 
mis-calibrat ion of the measurement devices, failure to perform 
prescribed data corrections, communications errors and simple 
recording errors. We present in this paper a quality control 
methodology to screen the data to identify those observations that 
appear to be erroneous in a statistical sense. It is then up to a 
meteorological analyst or a decision making algorithm to examine 
the tagged observations and decide on a course of action - delete 
the miscreant observations, modify the data if appropriate or 
include the observations, as is, for NOGAPS analysis. The 
procedures can also be used to monitor the reliability of the 
observing stations and to alert them of any mispract ices . 

Data Quality Tests: The methodology proposed here is an 
adaptation of a statistical test for outliers (Dixon, 1950 [1]). 
The first step in the application of this test is to partition the 
monitoring stations into small "contiguous groups". The recorded 
observations within any single group are expected to be "similar" 
because the weather patterns at these monitoring sites are 
comparable due to the proximity of the stations to each other. 
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One approach towards identifying such a group of stations is the 
following. Select the latitude and longitude for a target 
location on the surface of the earth (this could be the 

coordinates of a reporting station of interest) and choose an 

approximately square region by specifying a range for the 
latitude, say 5° above and below the target, and within a similar 
distance in the longitudinal direction. The observational 
increments (differences between the reported measurements and the 
forecast values) for all of the stations within the designated 
area are then subjected to a sequence of outlier tests. 

Assume that the number of reporting stations within the target 
area is N and the data from the most recent k days (say 30 days) 

is to be examined for quality. Let Xu , Xj. 2 , , Xi k be 

the observation increments for the i th station, i = 1, 2, . . . , 

N; let , M 2 , . . . , M n be -the means of the increments 

(averaged over the k days for each of the N stations) and Si, S 2 , 

. . . , S N the standard deviations (over the k days) for the N 

stations . 

The first test we propose is to compare the daily observational 
increments Xij , X 2 j , ... , X N j for the j th day, and check for 

"outliers". Let Y x < Y 2 < . . . < Y N be the data rearranged in 
increasing order. Compute 



r n _ (Yn ~ Yn_i) / (Y n - Y 2 ) ; 

if r n exceeds the tabulated value in Table A (extracted from 

Dixon, 1950) , corresponding to a chosen statistical significance 
level li (say .05), conclude that the observed datum for the 
station corresponding to Y N is an outlier. This implies that the 
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observation increment Y N is too large relative to the increments 
for the other stations in the region. That is, the observed 
geopotential height appears to be too high (large) when compared 
to the forecast . Now compute 



r'u = <Y 2 -Yi) / (Yn-x - Y x ) 

and conclude that Yx is an outlier if r'xx exceeds the tabulated 

value in Table A; if this happens, the implication is that the 
observed geopotential height is too low (small) compared to the 
forecast. These tests can be performed for each of the days for 
which data is available. The meteorological analyst or the 
decision making algorithm will then have to choose a course of 
action, based on some established decision rule, for each of the 
observations failing the test. The two tests were applied to 
actual data from three different areas of the globe viz., US, 
Europe, and China. In each area we identified a location by 

specifying its latitude and longitude and then selected all the 
stations within 4° of this location. The results of the daily 
tests are presented in Tables la,b,c. The I.D. numbers of the 
stations included for the test are listed in the first column. 
The entries (there are two entries in each of the 16 columns) in 
the table represent the number of days in October ' 90 on which the 
stations failed the outlier test at the low end (rxx' test) and the 
test at the high end (rxx test) using the data (geopotential height 

differences) for 16 pressure levels (1000, 850, 700, 500, 400, 

300, 250, 200, 150, 100, 70, 50, 30, 20, 15, 10 mb) . A dot in any 
position implies no failures; an asterisk indicates insufficient 
data to perform the tests. Such a table can be used in making 
quality judgments on the performance of the monitoring stations. 

Another test that can also be used to monitor the performance 
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quality of weather stations is to perform the two outlier tests on 
the averages Mx , M 2 , . . . , M N . The data for the tests are 

these averages rearranged in increasing order. Tables 2a,b,c 
present the same type of information as Tables la,b,c except that 
this time the entries are a 1 or a dot; a 1 signifies that the 
average for a station failed the outlier test while a dot would 
signify that the station passed the test. 

A third outlier test, that can be viewed as a test for 
consistency, is to check for an outlier in the station standard 
deviations Si , S 2 , . . . , S N . If Yx < Y 2 < . . . < Y N are 
these standard deviations in ordered form, compute 

r 10 = ( Yfj - Yfj_! ) / (Y n - Yx) 

and conclude that the standard deviation corresponding to Y N is 
statistically too large relative to the others, if rio exceeds the 

critical value in Table B (extracted from Dixon, 1950) at a 
significance level B (say .05) . Tables 3a,b,c are similar to 
Tables 1 and 2 using the station standard deviations, with the 
difference that only the high end outlier test is used. 

We have written a computer program, using Fortran, that will 
select (1) the reporting stations within a specified 
(approximately square) region around a specified location and (2) 
the appropriate data for these stations from the Naval 
Oceanographic and Atmospheric Research Laboratory (NOARL) files 
[2], and automatically performs all of the tests described above. 
Tables 1, 2 and 3, described above, are the actual outputs of this 
program. The program also generates additional information such 
as the number of times the daily outlier tests were performed 
(Table 4) and the number of days for which data was available 
(Table 5), during the period specified for the test. The numbers 
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in Table 4 will sometimes be smaller than the corresponding 
entries in Table 5; this happens whenever the number of data 
points available for the test is smaller than the minimum sample 
size of 4 . Table 6 is an example of another type of output of the 
program; for a selected station the table provides the daily 
test history, over the 16 pressure levels. A "PASS" indicates 
that the station passed all the daily tests; a "LFAIL" signifies 
that the observational increment is too low as compared to the 
numbers for the other stations in the region, on at least one day 
of the testing period. A "FAILH" denotes that the increment for 
at least one day was deemed to be high. Similar tables displaying 
the results of the means and standard deviations tests are also 
produced. In the "TOTAL" column a "PASS" means that the station 
passed all the tests on all the days and a "FAIL" implies the 
converse . 
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Extracted from: "Ratios of Extrema Values", W. J. Dixon, Annals of Mathematical 

Statistics. Vol. 22, No. 1, March, 1951. 



STATION / Lvl 1 2 3 4 5 



6 7 8 9 10 11 12 13 14 15 16 



72456 


♦ 


♦ 


, 




, 


• 


t 


2 




• 


1 




1 




1 


1 


3 


• 


t 




2 






# 


2 


# 




i 


i 


2 


* 


* 


72357 


• 


• 


2 




5 


• 


2 




1 


• 


1 


, 


2 


• 


1 


, 


1 


• 


1 




2 






# 


2 




1 




i 


• 


* 


* 


72349 


• 


1 


1 






















































• 


* 


* 


72340 
































1 


1 
























i 


• 


* 


* 


72260 


t 


2 




















1 


• 


1 


2 


1 


• 


t 




♦ 


, 






1 


• 


1 


• 




, 


2 


* 


* 


72247 


1 


, 


1 


1 




1 


• 


t 


• 


t 




, 


1 


♦ 




• 


1 


• 


• 


• 


1 




1 


• 


1 


, 


3 




5 




* 


* 


72239 


































* 


* 


* 


* 


* 


* 


* 


* 


* 


* 


* 


* 


* 


* 


* 


* 


72257 


* 


* 


* 


* 


* 


* 


* 


* 


* 


* 


* 


* 


* 


* 


* 


* 


* 


* 


* 


* 


* 


* 


* 


* 


* 


* 


* 


* 


* 


* 


* 


* 



Table la: Failures in daily outlier test (. =0, * = no test) 

U.S. stations near (35 ,-95 ) 

10/01/90 to 10/31/90 at 00Z 



STATION / Lvl 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 
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Table lb: Failures in daily outlier test (. =0, * = no test) 

European stations near (55°, 15°) 

10/01/90 to 10/31/90 at 00Z 



STATION / Lvl 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 
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Table lc: Failures in daily outlier test (. =0, * = no test) 

Chinese stations near (32°, 110°) 

10/01/90 to 10/31/90 at 12Z 
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STATION / Lvl 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 

72456 

72357 ....1.1 

72349 

72340 

72260 

72247 1 . 1 ... 

72239 ******************************** 

72257 ******************************** 

Table 2a: Failures in the means test (. =0, * = no test) 

U.S. stations near (35°, -95°) 

10/01/90 to 10/31/90 at 00Z 



STATION / Lvl 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 

26702 

26406 

12425 

12374 * * 

12330 **** 

12120 ' ****** 

10437 ********** 

10338 

10046 ********** 

10035 

9393 

9184 

6181 

2591 **** 

2527 **** 

Table 2b: Failures in the means test (. =0, * = no test) 

European stations near (55°, 15°) 

10/01/90 to 10/31/90 at 00Z 



STATION / Lvl 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 

57679 * * 

57516 * * . 1 1 1.1.1.. 

57494 

57461 ..** ** 

57447 . . * * 

57328 **** l , , ,****************** 

57178 ..** ** 

57127 

57083 * * 

57036 * * 

57245 **** ********** 

57633 ******************************** 

Table 2c: Failures in the means test (. =0, * = no test) 

Chinese stations near (32°, 110°) 

10/01/90 to 10/31/90 at 00Z 
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STATION / Lvl 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 

72456 

72357 

72349 

72340 

72260 

72247 1 ... 1 1 

72239 **************** 

72257 **************** 

Table 3a: Failures in the standard deviations test (. =0, * = no test) 

U.S. stations near (35°, -95°) 

10/01/90 to 10/31/90 at 00Z 



STATION / Lvl 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 

26702 

26406 

12425 1 

12374 * 

12330 .1 * * 

12120 ' *** 

10437 * * * * * 

10338 

10046 * * * * * 

10035 

9393 

9184 

6181 

2591 * * 

2527 * * 

Table 3b: Failures in the standard deviations test (. =0, * = no test) 

European stations near (55°, 15°) 

10/01/90 to 10/31/90 at 00Z 



STATION / Lvl 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 

57679 * 

57516 * 

57494 

57461 * * 

57447 * 

57328 * * ********* 

57178 * * 

57127 

57083 * 

57036 * 

57245 ** * * * * * 

57633 **************** 

Table 3c: Failures in the standard deviations test (. =0, * = no test) 

Chinese stations near (32°, 110°) 

10/01/90 to 10/31/90 at 00Z 
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0 


0 


0 


0 


0 




Table 4: 


Number of 


days 


; daily outlier 


test 


was 


performed 










U.S. 


stations 


near 


■ (35 


o 

' >“ 


■95°) 




















10/01/90 


to 


10/31/90 


at 


ooz 













STATION / 


Lvl 1 


2 


3 


4 


5 


6 


7 


8 


9 


10 


11 


12 


13 


14 


15 


72456 


31 


30 


31 


31 


31 


31 


31 


31 


31 


31 


30 


29 


28 


25 


24 


72357 


31 


30 


31 


31 


31 


30 


31 


31 


31 


30 


30 


30 


30 


28 


25 


72349 


31 


30 


31 


31 


30 


31 


30 


30 


30 


30 


30 


29 


28 


28 


23 


72340 


31 


29 


31 


31 


30 


30 


30 


30 


30 


30 


29 


28 


28 


27 


24 


72260 


30 


29 


30 


30 


30 


30 


30 


29 


29 


29 


29 


29 


29 


28 


27 


72247 


31 


28 


31 


31 


31 


31 


31 


'31 


31 


31 


29 


29 


29 


27 


23 


72239 


1 


1 


1 


1 


1 


1 


1 


1 


0 


0 


0 


0 


0 


0 


0 


72257 


0 


0 


0 


0 


0 


0 


0 


0 


0 


0 


0 


0 


0 


0 


0 



Table 5: Number of daily reports 

U.S. stations near (35°, -95°) 
in /m /on IP/'ll/Qn o + 007 

x w/ u x/ wJ iw/ J 1/ v»w ixi, Lt 



STATION 72456: 


LAT , LON 


= 39, -95 










LEVEL 


#DAYS 


DAILY TESTS 


MEANS TESTS 


STD TESTS 


TOTAL 


1 


31 


PASS 


PASS 


PASS 


PASS 


PASS 


PASS 


2 


30 


PASS 


PASS 


PASS 


PASS 


PASS 


PASS 


3 


31 


PASS 


PASS 


PASS 


PASS 


PASS 


PASS 


4 


31 


PASS 


*FAILH 


PASS 


PASS 


PASS 


*FAIL* 


5 


31 


PASS 


PASS 


PASS 


PASS 


PASS 


PASS 


6 


31 


LFAIL* 


PASS 


PASS 


PASS 


PASS 


*FAIL* 


7 


31 


LFAIL* 


PASS 


PASS 


PASS 


PASS 


♦FAIL* 


8 


31 


LFAIL* 


*FAI LH 


PASS 


PASS 


PASS 


*FAI L* 


9 


31 


LFAIL* 


PASS 


PASS 


PASS 


PASS 


♦FAIL* 


10 


31 


PASS 


PASS 


PASS 


PASS 


PASS 


PASS 


11 


30 


LFAIL* 


PASS 


PASS 


PASS 


PASS 


*FAIL* 


12 


29 


PASS 


PASS 


PASS 


PASS 


PASS 


PASS 


13 


28 


LFAIL* 


PASS 


PASS 


PASS 


PASS 


*FAIL* 


14 


25 


PASS 


*FAI LH 


PASS 


PASS 


PASS 


*FAIL* 


15 


24 


LFAIL* 


*FAI LH 


PASS 


PASS 


PASS 


♦FAIL* 


16 


11 


NOTEST 


NOTEST 


PASS 


PASS 


PASS 


PASS 






Table 6 


: Summary report for 


station 


72456 










10/01/90 to 


10/31/90 


at OOZ 







16 

0 

0 

0 

0 

0 

0 

0 

0 

16 

11 

4 

6 

2 

3 

6 

0 

0 



10 



non o o o 



Appendix: Program listing 



C THIS PROGRAM IS FOR CHECKING THE STANDARD FORECAST DEVIATION 

C FILES FOR OUTLIERS IN: (1) DAILY READINGS, (2) MEANS, AND 

C (3) STANDARD DEVIATIONS OF STATIONS WITHIN A TARGET DISTANCE 

C OF A GIVEN LATITUDE AND LONGITUDE. 

C 

C WRITTEN BY RICHARD FRANKE 

C DEPARTMENT OF MATHEMATICS 

C NAVAL POSTGRADUATE SCHOOL 

C MONTEREY, CA 93943 

C 0083P@CC.NPS.NAVY.MIL 

C 408/646-2758 

C 

C DATE OF BEGINNING: 6/21/91 

C DATES OF CHANGES: 6/21/91 

C MODIFIED TO USE SYMBOLS FOR OUTPUT: 6/21/91 

C MINOR PRETTYING UP: 6/24/91 

C MODIFIED TO NOTE NUMBER OF TESTS AND SUMMARY BY STATION: 6/27/911 

PAR.AMETER (MX = 152 ,NS=15 ) 

REAL FE(MX , NS , 16 ) , FER (MX , NS ) , XMN ( NS ) , STD ( NS) , FED ( NS ) 

INTEGER YRMD( MX , NS ) , NSTAT( NS ) , KNT (NS) ,NSD(NS) , 

1 NDLF ( NS , 16 ) , NDHF( NS, 16 ) , NSDF( NS ,16), NMNLF( NS , 16 ) , NMNHF ( NS, 16 ) , 

2 NSS(NS) ,YMD(MX,NS) ,NDG(NS, 16) ,NDDT(NS,16) ,LATS(NS) , LONS (NS) 
CHARACTER* 15 IUNIT.OUNIT 

CHARACTERS FLAGS (-1:1), FLAGSL (-1:1), FLAGSH (-1:1) 

CHARACTER* 1 SY( - 1 : 36 ) 

CHARACTER*38 SYY 
EQUIVALENCE (SY.SYY) 

DATA SYY/’ * . 1 234567 89ABCDEFGH I JKLMNOPQRSTUVWXYZ& ’ / 

DATA FLAGS/ ’NOTEST’ , ’ PASS \’*FAIL*’/, 

1 FLAGSL/’ NOTEST’ , ’ PASS ’,’LFAIL*’/, 

2 FLAGSH/ ’ NOTEST’ , ’ PASS ’,’*FAILH’/ 

GET DATA INPUT FILE NAME AND OUTPUT FILE NAME 
WRITE( * , 9 ) 

9 FORMAT! ’ NAME THE INPUT FILE’) 

READ ( * , 10 ) IUNI T 
10 FORMAT (A1 5) 

OPEN ( 4 , BLANK= ’ ZERO ’ , FI LE=I UNI T , MODE= ’ READ ’ , STATUS= ’ OLD ’ ) 

WRI TE( * , 21 ) 

21 FORMAT! ’ NAME THE OUTPUT FILE’) 

READ!*, 10)OUNIT 

OPEN ( 8 , FI LE=OUNIT , MODE= ’WRITE ’ , STATUS= ’ NEW ’ ) 

GET THE INPUT 
100 WRITE!*, 1) 

1 FORMAT!’ INPUT THE TARGET LATITUDE, LONGITUDE AND THE RANGE’,/ 

1 ’ - FREE FIELD - RANGE = 0 STOPS EXECUTION’,/) 

READ! *,*)LAT, LON, LLR 
I F ( LLR. LE. 0 )STOP 
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STRETCH OUT IN THE LONGITUDINAL DIRECTION 

LTR=LLR/COS ( FLOAT ( LAT ) /57 . 3 ) 

START INPUT UNIT 4 AT THE BEGINNING 

REWIND 4 
MS = 0 



READ A HEADER 

200 READ (4,2, END=300 ) NST , I LAT , I LON , I TM , NREC 

2 FORMAT (BZ,3X,I7,4X,2I5,8X,I2,I5) 

I F ( ABS ( I LAT-LAT ) . LE . LLR .AND. ABS ( I LON-LON ) . LE . LTR ) THEN 

IT’S IN THE RECTANGLE, GET THE DATA 

MS=MS + 1 
NSTAT(MS)=NST 
LATS(MS )=ILAT 
LONS(MS)=ILON 

READ ( 4 , 3 ) ( { FE ( I , MS ,J) ,J=1,16) , YRMD( I ,MS) , 1=1 , NREC ) 

3 FORMAT( BZ , 16F4 .0,18) 

ELSE 

OR SKIP IT 

READ (4,4)(JUNK,I = 1, NREC ) 

ENDIF 

4 FORMAT (II) 

GO BACK FOR MORE IF WE DON’T HAVE THE MAX 

IF(MS.LT.NS) GO TO 200 
PRINT 7 

7 FORMAT (’ STOPPING BEFORE END’) 

300 CONTINUE 
IF(MS. LT. 4 ) THEN 

WRITE(*,5)MS 

NOT ENOUGH DATA, START OVER 

5 FORMAT ( ’ ONLY’, 1 2,’ STATIONS FOUND - START OVER’) 

GO TO 100 

ELSE 

WRITE( *, 6 )MS, ( NSTAT ( N ) , N=1 ,MS) 

6 FORMAT (15,’ STATIONS FOUND: ’,/( 1 10 ) ) 

ENDIF 

GET THE DATE RANGE 

301 WRITE( * , 8 ) 
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8 FORMAT ( ’ INPUT THE BEGINNING AND ENDING DAYS (> = 1, < = 152 ) ’ / 

1 ’ BEGINNING = ZERO GOES TO NEW LATITUDE/LONGITUDE’/ 

2 ’ BEGINNING = NEGATIVE ENDS RUN’) 

READ ( * , * ) I DB, IDE 
IF(IDB.LT. 0)STOP 
I F( IDB. EQ. 0 )GO TO 100 

OK, LET’S DO THE JOB 

ZERO OUT FLAG FIELDS 

DO 304 L=1 , 16 
DO 303 N=1 , NS 
NDLF(N,L) = 0 
NDHF( N, L ) = 0 
NSDF (N , L ) = -1 
NMNLF(N, L) = -1 
NMNHF( N , L ) = -1 
NDG( N , L ) = 0 
NDDT( N , L ) = 0 

303 CONTINUE 

304 CONTINUE 

LOOP OVER ALL LEVELS 

DO 400 L=1 , 16 
KS=0 

KS COUNTS STATIONS WITH A HISTORY OF 2 OR MORE DAYS 

DO 305 N=1 ,MS 
KNT(N) = 0 
XMN(N) = 0. 

KNT COUNTS HOW MANY DAYS 

305 CONTINUE 
DO 310 I=IDB, IDE 

ND=0 

DO 308 N=1 , MS 
C 

C ND COUNTS THE NUMBER OF STATIONS KITH GOOD DATA TODAY 

C FED SAVES IT FOR THE TEST 

C NSD TELLS WHICH STATIONS HAVE GOOD DATA 

C FER SAVES IT ALL FOR MEAN AND STD CALCS 

C YRMD SAVES THE DATE INFORMATION (NOT USED PRESENTLY, 6/21 

C 

FED( N)=0. 

IF(FE(I ,N,L). LT .999) THEN 
NDG ( N , L ) =NDG ( N , L ) + 1 
KNT(N) = KNT(N) + 1 
ND=ND+1 

FED(ND)=FE( I ,N,L) 

NSD( ND)=N 
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FER( KNT ( N ) , N ) = FE(I,N,L) 

YMD(KNT( N ) , N ) = YRMD( I , N ) 

XMN(N)=XMN(N)+FER(KNT(N),N) 

ENDIF 

308 CONTINUE 

TEST CURRENT DAY’S DATA 

CALL R1 1 P05 ( FED , ND , NB , I B , NT , I T , I ER ) 

IF(IER.LE.O) THEN 

IF(NB.GT.O)NDLF(NSD(IB) ,L) = NDLF ( NSD( IB) , L ) + 1 
I F( NT. GT. 0)NDHF(NSD( IT ) , L ) = NDHF(NSD( IT) ,L) + 1 
DO 309 N=1 , ND 

NDDT(NSD(N) ,L) = NDDT(NSD(N) ,L) + 1 

309 CONTINUE 
ENDIF 

310 CONTINUE 

DO 315 N=1,MS 

IF(NDDT(N,L) .EQ.O) THEN 
NDLF(N,L) = -1 
NDHF( N, L) = -1 
ENDIF 

315 CONTINUE 

DO 330 N=1 , MS 

IF(KNT(N).GT. DTHEN 

CALCULATE MEAN AND ST. DEV IF MORE THAN ONE DAY OF DATA 

KS=KS+1 

NSS(KS) = N 

XMN ( KS ) =XMN ( N ) /KNT ( N ) 

S= 0. 

DO 320 1=1 , KNT ( N ) 

S=S+ ( FER ( I , N ) -XMN ( KS ) ) **2 
320 CONTINUE 

STD(KS) = SQRT ( S/ ( KNT ( N ) - 1 ) ) 

ENDIF 

330 CONTINUE 

RUN TEST ON MEAN VALUES FOR THE STATIONS 

CALL R11P05 ( XMN , KS , NB, IB, NT , IT , IER ) 

IF( IER. EQ.O) THEN 
DO 340 N=1 , KS 
NMNLF(NSS(N),L) = 0 
NMNHF(NSS(N ) ,L) = 0 
340 CONTINUE 

IF(NB.NE.O) NMNLF(NSS(IB) ,L)=1 
IF(NT.NE.O) NMNHF(NSS(IT),L)=1 
ENDIF 

RUN TEST ON STANDARD DEVIATIONS OF THE STATIONS 
CALL R10P05(STD,KS,NT,IT,IER) 
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I F( I ER. EQ. 0 ) THEN 
DO 360 N=1 ,KS 

NSDF( NSS( N ) , L ) = 0 
360 CONTINUE 

IF(NT.GT.O)NSDF(NSS(IT),L) = 1 
ENDIF 

400 CONTINUE 
MDATB=0 
MDATE=0 
DO 410 N=1 ,MS 

I F ( MDATB . EQ . 0 ) MDATB= YRMD ( I DB , N ) 

I F ( MDATE . EQ . 0 ) MDATE=YRMD ( I DE , N ) 

410 CONTINUE 

I F ( MDATB . EQ . 0 )MDATB=MDATE- ( I DE- I DB ) 

I F ( MDATB . LT . 0 ) MDATB= I DB 
I F ( MDATE . EQ . 0 ) MDATE=MDATB+ ( I DE- I DB ) 

OUTPUT RESULTS FOR VARIOUS TESTS 

WR I TE ( 8 , 1 1 ) MDATB , MDATE , I TM , L AT , LON , LLR , LTR 

11 FORMAT! ’ 1 RESULTS OF THE OUTLIER TESTS FROM’, 

1 18, ’ TO’ ,18, ’ AT’ ,13, ’ Z ’// 

2 ’ NUMBER OF FAILURES BY LEVEL AND STATION: .=0, *=NOTEST’ 

3 //’ TARGET LATITUDE/LONGITUDE IS’ ,215, ’, RANGE: ’ ,12, ’X’ ,12 

4 //’ RESULTS OF THE DAILY OUTLIER TEST (LOW FAIL/HIGH FAIL)’) 
WR I TE ( 8 , 1 2 ) ( ( L ) , L= 1 , 1 6 ) 

12 FORMAT! / ’ STATION / Lvl’,1614/) 

WRITE (8, 13) (NSTAT(N) , 

1 ( SY ( NDLF ( N , L ) ) , SY( NDHF( N , L ) ) , L=1 , 16 ) , N=1 ,MS ) 

13 FORMAT! 18, 5X , 32A2 ) 

WRITE! 8 , 20 ) IDE-IDB+1 

20 FORMAT!//’ NUMBER OF DAYS DAILY OUTLIER TEST WAS PERFORMED’ 

1 ’ OUT OF’ , 14) 

WRITE (8, 12)((L),L= 1,16) 

WRITE! 8, 18) (NSTAT(N) , (NDDT(N,L) , L=1 , 16 ) ,N=1,MS) 

18 FORMAT! 18, 4X, 1614) 

WRITE! 8,14) 

14 FORMAT!//’ RESULTS OF THE MEANS OUTLIER TEST (LOW FAIL/HIGH’ 

1 ’FAIL)’) 

WRITE! 8 , 12 ) ( ( L ) , L=1 , 16 ) 

WRITE! 8, 13) (NSTAT(N) , 

1 ( SY(NMNLF(N , L) ) , SY( NMNHF! N , L ) ) , L= 1,16) ,N=1,MS) 

WRITE (8 , 15 ) 

15 FORMAT!//’ RESULTS OF THE STANDARD DEVIATION OUTLIER TEST’) 
WRITE! 8, 12 ) ( ( L) , L=1 , 16 ) 

WRITE! 8 , 16 ) ( NSTAT ( N ) , ( SY ( NSDF ( N , L ) ) , L=1 ,16), N= 1 , MS ) 

16 FORMAT! 1 8 , 4X , 16A4 ) 

WR I TE ( 8 , 1 7 ) I DE- I DB+ 1 

17 FORMAT!//’ NUMBER OF DAYS REPORTING OUT OF’ ,14) 

WRITE (8, 12)((L),L= 1,16) 

WRITE! 8, 18) ( NSTAT! N) , (NDG(N, L) ,L=1 , 16 ) ,N=1,MS) 

OUTPUT SUMMARY BY STATION 
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c 

c 

c 

c 

c 

c 

c 



WRITE( 8 ,22 ) 

22 FORMAT(//’ ************’, 

2 //’ SUMMARY BY STATION FOLLOWS’/) 

DO 500 N=1,MS 

WRITE(8,23)NSTAT(N),LATS(N) ,LONS(N)- 

23 FORMAT (//’ 

1 ’ 

2 /’ STATION’ ,16, ’ : LAT,LON = ’,14, ’,’,14/ 

3 /’ LEVEL #DAYS DAILY TESTS MEANS TESTS 

4 ’ STD TESTS TOTAL’/) 

DO 420 L=1 , 16 

SET CODES FOR TESTS: -1 IS NOTEST, 0 IS PASS ALL, 1 IS FAIL 

ONE OR MORE 

FOR THE TOTAL TEST, NOTEST MEANS NO TESTS PERFORMED, 

PASS MEANS PASS ALL TESTS PERFORMED (IF ANY) 
FAIL MEANS FAILED AT LEAST ONE 

NFL=MIN(NDLF(N,L) ,1) 

NFH=MI N ( NDHF ( N , L ) , 1 ) 

NMFL=MI N ( NMNLF ( N , L ) , 1 ) 

NMFH=MI N ( NMNHF ( N , L ) , 1 ) 

NSF=MIN(NSDF(N,L) ,1) 

NTF=MAX ( NFL , NFH , NMFL , NMFH , NSF ) 

WRITE( 8 , 24 )L , NDG( N, L ) , FLAGSL( NFL ) , FLAGSH( NFH ) , 

1 FLAGSL ( NMFL ) , FLAGSH ( NMFH ) , FLAGS ( NSF ) , FLAGS (NTF) 

420 CONTINUE 

24 FORMAT ( 21 6 , 4X , A6 , 3X , A6 , 5X , A6 , 3X , A6 , 5X , A6 , 5X , A6 ) 

500 CONTINUE 

WRITE( 8,19) 

19 FORMAT!//' ******************’, 

1 ’ ********************************************************* ’/////) 
GO TO 301 
END 



SUBROUTI NE R10P05 ( X , N , NOUTT , IT , I ER ) 

PARAMETER (MAX=30) 

DIMENSION X ( N ) , SX ( MAX ) , I SX ( MAX ) 

DIMENSION CV( 30 ) 

C 

C THIS SUBROUTINE DOES THE RIO TEST FOR OUTLIERS AT THE .05 CRITICAL 
C TEST LEVEL. 

C 

C THE CV VALUES ARE THE CRITICAL VALUES FROM TABLE V. 

C 

DATA CV/2*0, .941, .765,. 642, .560, .507,. 468, .437, .412, .392, .376, 

1 .361, .349, .338, .329, .320, .313, .306, .300, .295, .290, .285, .281, 

2 .277, .273, .269, .266, .263, .260/ 

C 

C INPUT VALUES ARE: X - THE DATA TO BE TESTED 

C N - NUMBER OF DATA POINTS 

C 

C OUTPUT VALUES ARE: NOUTT - FLAG FOR LARGEST VALUE 
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c 

C IB - 

C IER 

C 

C 

C 

C 

C 

C CHECK FOR SUITABLE N 
C 

IER = 0 
NOUTT = 0 

I F ( N . LT . 3 . OR . N . GT . MAX ) THEN 
IER = 1 
RETURN 
END I F 

DO 100 1=1, N 
SX ( I ) = XU) 

ISX(I ) = I 
100 CONTINUE 



0 IF OK, 1 IF NOT 
LOCATION IN X ARRAY 
ERROR RETURN INDICATOR 

0 IF ALL OK 

1 IF N OUT OF RANGE 4 TO MAX, INCLUSIV 

NO TEST PERFORMED 



SORT INTO INCREASING ORDER 

CALL SHSORT ( SX , ISX , N ) 

RT = 0. 

I F( SX ( N)-SX ( 1 ) .NE. 0.) RT = ( SX ( N ) - SX(N-l) )/(SX(N) - SX(1)) 
IF(RT.GT.CV(N) )THEN 
NOUTT = 1 
IT = ISX(N) 

ENDIF 

RETURN 

END 



C 

C 

C 

C 

C 

C 



C 

C 

C 

C 

C 

C 

C 

C 

c 



SUBROUTINE R1 1P05 ( X , N , NOUTB, I B , NOUTT , I T , I ER ) 

DIMENSION X ( N ) , SX ( 30 ) , ISX ( 30 ) 

DIMENSION CV ( 30 ) 

THIS SUBROUTINE DOES THE Rll TEST FOR OUTLIERS AT THE .05 CRITICAL 
TEST LEVEL. 

THE CV VALUES ARE THE CRITICAL VALUES FROM TABLE V. 

DATA CV/3*0, .955, .807, .689, .610, .554, .512, .477, .450, .428, .410, 

1 .395, .381, .369, .359, .349, .341, .334, .327, .320, .314, .309, .304, 

1 .299, .295, .291, .287, .283/ 

INPUT VALUES ARE: X - THE DATA TO BE TESTED 

N - NUMBER OF DATA POINTS 

OUTPUT VALUES ARE: NOUTB - FLAG FOR SMALLEST VALUE 

0 IF OK, 1 IF NOT 
IB - LOCATION IN X ARRAY 
NOUTT - FLAG FOR LARGEST VALUE 
0 IF OK, 1 IF NOT 
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C IT - LOCATION IN X ARRAY 

C IER ERROR RETURN INDICATOR 

C 0 IF ALL OK 

C 1 IF N OUT OF RANGE 5 TO 30, INCLUSIVE 

C 

C 

C CHECK FOR SUITABLE N 
C 

IER = 0 

I F ( N . LT . 4 . OR . N . GT . 30 ) THEN 
IER = 1 
RETURN 
ENDIF 

DO 100 1=1, N 
SX(I) = X(I) 

ISX(I) = I 
100 CONTINUE 

SORT INTO INCREASING ORDER 

CALL SHSORT ( SX , I SX , N ) 

RB = 0. 

IF(SX(N-1 )-SX( 1 ) .NE. 0.) RB = ( SX ( 2 ) - SX( 1 ) )/(SX(N-l ) - SX ( 1 ) ) 
NOUTB = 0 

IF(RB.GT.CV(N) )THEN 
NOUTB = 1 
IB = ISX(l) 

ENDIF 
RT = 0. 

IF(SX(N)-SX(2) .NE. 0.) RT = (SX(N) - SX(N-l) )/(SX(N) - SX( 2 ) ) 
NOUTT = 0 

I F ( RT . GT . CV ( N ) )THEN 
NOUTT = 1 
IT = ISX(N) 

ENDIF 
RETURN 
END 

C Included for completeness oniy - from NONIMSL library at NPS 



C 

c 

C A. IDENTIFICATION: 

C TITLE: NUMERICAL SORT 

C ID: Ml-NPG-SHSORT (F-IV) 

C PROGRAMMER: R. BRUNELL 

C DATE: MARCH 1968 

C MODIFIED: DEC. 1973 BY L. NOLAN 

C 

C B. PURPOSE: 

C TO SORT, IN ASCENDING ORDER, AN ARRAY OF SINGLE PRECISION REAL 

C NUMBERS BY THE METHOD OF SHELL, AND TO PRODUCE AN ARRAY OF INDEXES 

C SO USER CAN RE-ORDER OTHER CORRESPONDING INFORMATION ACCORDING TO 
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c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 



ASCENDING VALUES OF "A". 

C. USAGE: 

1. CALLING STATEMENT: 

CALL SHSORT ( A , KEY , N ) 

2. ARGUMENTS: 

A - ARRAY OF NUMBERS TO BE SORTED. THIS ARRAY IS SORTED 
(RE-ORDERED) BY "SHSORT". 

KEY - ARRAY, DIMENSIONED AT LEAST N IN CALLING PROGRAM, TO BE 
FILLED BY USER WITH INTEGERS FROM 1 TO N. AFTER EXIT 
FROM SHSORT, KEY ( 1 ) WILL CONTAIN THE ORIGINAL INDEX OF 
THE SMALLEST ELEMENT OF "A"; KEY ( 2 ) WILL CONTAIN THE 
ORIGINAL INDEX OF THE NEXT-TO-SMALLEST ELEMENT OF "A"; 
ETC. KEY ( N ) WILL CONTAIN THE ORIGINAL INDEX OF THE 
LARGEST ELEMENT OF "A". 

N - NUMBER OF MEMBERS IN ARRAYS "A" AND "KEY". 

D. REFERENCES: 

1. "ALGORITHM 201, SHELLSORT" , BOOTHROYD, J. , "COMMUNICATIONS OF 

ACM", VOL 6, NO 8, AUGUST 1963, P.445. 

2. "CERTIFICATION OF ALGORITHM 201", BATTY, M. A., "COMMUNICATIONS 

OF ACM", VOL 7, NO 6, JUNE 1964, P.349. 

SUBROUTINE SHSORT ( A , KEY , N ) 

DIMENSION A ( N ) , KEY ( N ) 

Ml = l 

6 M1=M1*2 

IF (Ml .LE. N) GO TO 6 
Ml=Ml/2-l 
MM=MAX0(Ml/2, 1) 

GO TO 21 

20 MM=MM/2 

IF (MM .LE. 0) GO TO 100 

21 K=N-MM 

22 DO 1 J=1 , K 
1 1 =J 

11 IM=I I +MM 

IF ( A ( I M ) .GE. A ( 1 1 ) ) GO TO 1 
TEMP=A( 1 1 ) 

IT=KEY( 1 1 ) 

A ( 1 1 )=A( IM) 

KEY( 1 1 )=KEY( IM) 

A(IM)=TEMP 
KEY ( I M ) = I T 
1 1 = 1 1-MM 

IF (II . GT . 0 ) GO TO 11 
1 CONTINUE 
GO TO 20 
100 RETURN 
END 
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